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Abstract 


In the current study, a new numerical algorithm is presented to solve a class 
of nonlinear fractional integral-differential equations with weakly singular 
kernels. Cubic hat functions (CHFs) and their properties are introduced for 
the first time. A new fractional-order operational matrix of integration via 
CHFs is presented. Utilizing the operational matrices of CHFs, the main 
problem is transformed into a number of trivariate polynomial equations. 
Error analysis and the convergence of the proposed method are evaluated, 
and the convergence rate is addressed. Ultimately, three examples are 
provided to illustrate the precision and capabilities of this algorithm. The 
numerical results are presented in some tables and figures. 
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1 Introduction 


The mathematical modeling of many phenomena in various branches of sci- 
ence leads to nonlinear integral-differential equations. Fractional calculus is 
applied extensively by many scientists in the mathematical modeling and 
control of numerous dynamic systems [30, 31]. This class of equations arises 
in the field of signal processing [21], waves and brain modeling [18, 26], ra- 
diative equilibrium [17], and so on. Commonly, it is impractical to obtain an 
analytical solution to integral-differential equations. As a result, the improve- 
ment of some numerical methods and the introduction of new high-accuracy 
numerical algorithms is very important to obtain approximate solutions. So 
various numerical methods have been developed to solve these types of equa- 
tions by many researchers. Some of the prominent methods are modified 
differential transform [15, 20], Adomian decomposition, Homotopy analysis 
(9, 12], Galerkin [10, 28], collocation [16, 25], product integration [1], Euler 
wavelets [7], haar wavelets [4], Legendre wavelets [29], Chebyshev wavelets 
[13], Hermite cubic splines [27], Hat functions [11], Taylor series [8], and so 
forth. The nonlinear fractional integral-differential equation with a weakly 
singular kernel appears in the following form: 


C D2u(t) = g(t) + acer f( (t—s) °u™(s)ds, a>0, 0<B<1, 


ates i=0,1,...,Jfa]-1, meéeN, tel(t), 
(1) 


where u(t) is an unknown function to be determined, » is an appropriate 
parameter, g(t) and p(t) are known continuous functions on I(T) := [0,7], 
and ( D@ is the Caputo fractional differential operator of order a. Some nu- 
merical methods convert such an integral-differential equation into a system 
of algebraic equations that can be easily solved. 
Wang and Zhu [32] applied the second kind of Chebyshev wavelets method 
to give approximate solutions for the fractional integral-differential equations 
with a weakly singular kernel. Nemati and Lima [22] applied a numerical 
method based on modified hat functions (MHFs) for solving the problem (1). 
Xie et al. [33] used the Haar wavelets to solve a coupled system of fractional- 
order integral-differential equations. Riahi Beni [29] proposed a novel tech- 
nique for nonlinear fractional Volterra—Fredholm integro-differential equa- 
tions. Also, a numerical solution for a fractional integro-differential equation 
via a method based on the Gegenbauer wavelets was suggested by Ozaltun, 
Konuralp, and Giimgiim [23]. In this paper, we introduce a high-precision 
numerical algorithm for the problem (1) in terms of Cubic hat-functions 
(CHFs). 

The present work discusses some of the properties of Riemann—Liouville 
integral operators to solve the nonlinear fractional integral-differential equa- 
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tions. Applying the operational matrix method, the principal problem will 
be reduced to solving several nonlinear trivariate polynomial equations. In 
Section 2, some basic definitions and characteristics of fractional calculus 
are presented. Section 3 is devoted to introducing the operational matrix of 
CHF’ basis. The fourth section studies the absolute error of approximation 
of a function by a truncated series of CHFs. The fifth section presents a nu- 
merical method for Problem (1). The convergence analysis of the proposed 
scheme is discussed in Section 6. To show the validity and accuracy of the 
utilized approach, three numerical examples are provided in Section 7, and 
the paper ends in Section 8, with a conclusion and discussion. 


2 Basic concepts and definitions 


In this section, some definitions and properties, which have been used in this 
manuscript, are explained. In this research, the Riemann-—Liouville integral 
operator of the ath order (J?) and the Caputo fractional differential operator 
of order a (§.D#) will be used. They are well addressed in [24]. 


Definition 1. Suppose that ae R,n-1l<a<n,nEN, and let u(t) be 
a continuous function defined on [0,1]. The Caputo fractional derivative of 
order a > 0 is defined as follows: 


Tana) is (t — ye DF alr) dr, n-l<a<n, 


Cpa 
DFu(t) = 
Ot ( ) ul (t), a=n, 


(2) 


wherein 


re) = | es i 
0 


Definition 2. Assume that a > 0 and that u(t) is a continuous function 
defined on the closed interval [0, 1]. The Riemann-—Liouville integral operator 
of order a is defined as follows: 


Ieu(t) = Ta) | (t — 7)°"*u(r)dr. (3) 


The Riemann-Liouville integral operator and the Caputo fractional deriva- 
tive operator satisfy the following properties [24]: 


Fe (UPu(t)) = fut) = 17 tP ult), a, 8 > 0, 
n-1 


: ti 
Ip (§ Deu(t)) = u(t) — So uO) 5, p—-leesn, FS0. 


i=0 
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2.1 Definition of CHFs 


First, let us state a history of Hat functions and then some definitions and 
properties of CHFs. In 2011, Babolian and Mordad [6] described gener- 
alized Hat functions to solve systems of Fredholm and Volterra integral 
equations. In 2016, Mirzaee and Hadadiyan [19] introduced MHFs to solve 
Volterra—-Fredholm integral equations. In this paper, we improve the hat 
functions method and use the method for solving linear and nonlinear frac- 
tional integral-differential equations with weakly singular kernels. CHFs are 
defined on the closed interval [0,7] and have a hat-like shape. The inter- 
val is divided into n subintervals, with equal lengths h, where h = fo and 
n=3K, K EN. 

CHF’ are defined as follows: 


L(t—h)(t—2h)(t-— 3h), 0<¢< 3h, 
do(t) = anole — hy X (5) 
0, otherwise. 
For i = 3v — 2, v = 1,2,...,n/3, 
epi “| gas (t— (6— DAY(E~ (E+ DANC— E+IA), GF“ Yh SES G42, 6) 
0, otherwise. 
For i = 3v—1, vy =1,2,...,n/3, 
1 3 : a i A 
$i(t) -| ays (t— (4 2)A)(E— (4 — Dh)(E- (A+ Ih), ( a)hsts( +1)h, (7) 
0, otherwise. 
When i = 3v, vy = 1,2,...,(n — 3)/3, 
sis (t — (i — 3)h)(¢— (é-2)h)(t-—(@- Dh), (G3) SES IA, 
di(t) = gas (t— (i+ U)A)(E- (+ 2)h)E-— (6+. 3)h), th<t<(i+3)h, (8) 
0, otherwise, 
and 
‘r= { gas (t — (n — 3)h)(t— (n — 2)h)(t—(n—1)h), (n— ae <t<nh, (9) 
0, otherwise. 
A function u(t) can be expressed in terms of CHF’ as follows: 
u(t) © Un(t) = — ai@i(t) = ATH(t) = (4) A, (10) 
i=0 
so that : 
and 
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A= [ao,a1,..., An], (12) 


wherein a; = u(th), i = 0,...,”, are unknown coefficients of the CHFs. Fig- 
ure 1 shows the CHFs plotted on the interval [0, 1] for n = 6 using MATLAB 
package. 
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2.1.1 Properties of CHFs 


Using the CHFs definition, the following properties can be obtained: 
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Figure 1: Plots of the CHFs, up ton =6,T =1 


Salt) =1,  di(jh) = ‘ i (13) 
i=0 ’ : 


Multiplying both sides of this summation to ¢;(t) attains 


(>: asta) = 6;(t). (14) 
i=0 


Thus, for t = jh, we have 


3? dy(5h)dilsh) = o4(h), 


[(ds(ih)bo( ih) + +++ + (by sh)by (sh) + ++ + (b5(db)bn GR))] = VR), 
[(0j (Gh) x 0) +--+ + (bj (GR) x O|(GR)) +--+ + (Gj (Gh) x 0)] = aaa 
15 
As a result, 
bj (Gh); (Gh) = $5 (9h). (16) 
Taking these properties, one has 
bi(Gh)o;(Gh) © 0, fea (17) 


Then, from the relations (17 ) and (11), it can be concluded that 


@(t) 7 (t) ~ diag [do(t), d1(t),---,n-1(t), On(#)) = diag (®(é)). (18) 


2.1.2 Nonlinear approximation of CHFs 


Using (18) and (10), u(t) , m = 1,2,..., can be calculated as follows: 


u?(t) ~ A’ ®(t)67(t)A = A” diag(&(t)) A = A” diag(A)®(t) 
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u3(t) ~ w2(t)u(t) = AP B(t)@7 (t)A = AT diag(®(t))A = AP diag(A)®(t) 


= Aj ®(t), As = [6 @4, ae elle 


n 


u(t) ~ S > ai 4i(t) = AT O(t), Am =(a™,a%,...,a™?. 


i=0 


3 Operational matrices of CHFs 


(19) 


In this part of the study, we achieve the fractional-order integral operational 


matrix using CHFs. 


3.1 Fractional order operational matrix of integration 


Let us state the following theorem. 


Theorem 1. Let &(t) be given by (11) and let a > 0. Then 


TP O(t) ~ Q*@(t), 


(20) 


where Q® is called the (n + 1) x (n + 1) operational matrix of fractional 


integration of order a and is defined as follows: 


0 pi p2 P3 Pa Ps Po °°: 
0 01 02 03 04 05 06 °°° 
0 ky Ko K3 Ka K5 Ke te: 
O pa He M3 Ma Hs He + 
00 0 0 01, 62 03 ::: 


Qo = 0 0 0 0 Ki K2 K3 °°° 
6ia+4) | 0 0 OO yuo fg 
Qi: f ff 4 
000000 0 
000000 0 
000000 0 
wherein 


Pn—-2 Pn—-1 Pn 
On-2 On-1 On 
Kn-2 Kn-1 Kn 
Un-2 Pn-1 Ln 
On—5 In—4 On-3 
Kn—-5 Kn—4 Kn—-3 
Un—5 Hn—4 Ln-3 


O71 02 03 
Kl K2 K3 
Mi M2 b3 


pr = 6k%(a + 3)(a + 2)(a + 1) — 11k°11 (a + 8)(a@ + 2) 


+ 12k? (a +3) — 6h", &=1,2,3, 


(21) 
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pr = 6k°(a + 3)(a + 2)(a + 1) — (11k — 2(¢ —3)°*) (a +3)(a +2) 
ket? + (k — 3°") (a +3) 


( 
ee aie a) ae 


6 
on =3 (6(K)° "(a 4 3)(a +2) — 10(k)°*?(a + 3) + 6(n)°**) , k=1,2,3, 


on = 9 (2(k)°t! — (k att) (a +3)(a +2) 


Ke = —3 (3(k)** (a + 8)(@ + 2) — 8(k)**2( + 3) + 6(K)°**) , k=1,2,3, 
ig ((aer i hae (a + 3)(a +2) 
+6 (40°? +5(k — ay) (a +3) 
— 18 ( nets — (k—3)°**), ee 
be = 2(k)°** (a + 3)(a + 2) — 6(k)°F? (a + 3) + 6(K)OTF, =k =1,2,3, 


ie (Gar =i Gee 7) (a + 3)(a + 2) — 6(k)*+?(a + 3) 
+6 ((*? Ak — a) ,  k=4,5,6, 
tp = 2 (Oe tiggea\ aie 6)°**) (a +3)(a +2) 
6 (ne? (i= 6°) (a +3) 
+6 (a? 2(k — 3)°*3 + (k 6)°**) » BS Tesch (22) 


Proof. First, for ¢;(t), i = 0,...,n, we have the definition of the Riemann— 
Liouville integral operator as follows: 


il : 4 
I? ¢,(t) = —~ t—7)* o;(r)dr. 23 
PO) = By f 7) Mealaer (23) 
We expand Jf'¢;(t), in terms of the cubic hat basis functions as follows: 
Te o(t) ~ Yo vs6;(t), 1=0,...,n, (24) 
j=0 


where the values of I?'¢;(t) at jth node point, (jh), represent the coefficients 
ij. Thus, we have 
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a er ee 7” 
Vij = ol. (Gh—T)” “di(r)dr, i,7 =0,1,...,n. (25) 
Using equations (5)—(9), we calculate the integral (25). For 7 = 0, by substi- 
tuting (5) in (25), we introduce the coefficient as follows: 
6jo+8 — 19je+2 


411j°+1(a + 3) 
—6j*%(a+3)(a4+2)(a+1), 


105 ~~ 6T (a +4) 6 Cas = ayer) 
6 Cie +G- aye?) (a +3) 
i aie meee ae (a + 3)(a +2) 
—6j%(a + 3)(a + 2)(a +1), 


For i = 3v — 2, v = 1,2,...,n/3, we obtain 


0, jgsi-l, 
6(j —i+.1)°*8 
—10(j —i+1)°*?(a 4+ 3) tS 7 S242, 
+6(j —i+1)°*1(a + 3)(a + 2) 
y= 6((j-é+ 1°"? - G—-4-2)°**) 
2T(a + 4) _9 (50 ox ij" 


+4(j -1-2)**) (a +3) - FAs. 
43 (2G -i+ cen 
—(j-4-1)***) (a +3)(0 +2) 


(27) 
For i = 3v—1, v= 1,2,...,n/3, replacing (7) into Eq. (25) yields 
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0, 7x 45—2, 
6(j —i+.2)°+8 
—8(j —i+2)°7(a +3) ; G— 1g St, 
+3(j —i+ 2)°*1(a + 3)(a + 2) 
he 6((j-t+2)7*9 — (Gj -4-)7**) 
~~ Ora +4) _9 (46 4a)? 
+5(j -i-1)**”) (o+3) |, 7>a+2. 


+3((j-i+1)°" 
SS i) (a +3)(a +2) 


(28) 
Now, we attain (25) for i = 3v, v = 1,2,...,n/3, 
0, j<i-3, 
6(j (fj —4+3)9+8 
Fe a : t-2<j <1, 
2(j —i+3)° (a + 3)(a +2) 
(G-i+3) ‘aaae — 25 - i)***) 
—6(j ena oe : : 
‘ oe 
a((j-i+3) ) a+ S79 S04 ) 
vai = WaT -11(j “rs (a +3)(a +2) 
(G G13)" 
~2j — 2 a aay 
-6 ((j-i+3)°? 
—(j-i-3)"?) (a +3) ee ees 
(j-—i4+3)0+? 
+2] —11(j-i)%*! | (a+3)(a 4 2) 
+(j —i-3)°*1 
(29) 


Consider 3v — 2 = 7 in (27), 3v — 1 = 7 in (28), and 3v = i in (29). Then 
apply 3v +k = 7 to all (26)-(29), v = 1,...,n/3 and k = 1,...,n. Some 
simple manipulations completes the proof. 


As a result of using (10) and (20), we can approximate the integral of a 
nonlinear function as follows: 
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Ieu™(t) ~ Ie (>: o)) ~ If (AT O(t)) ~ ATQ°H(t), m= 1,2,.... 
1=0 


(30) 
For instance, when a = 1 and n = 3, using the operational matrix (21), we 
get 


Examples Composite trapezoidal rule CHFs 

solutions with n = 6 solutions 
Example 1 :f sins cos’s ds 0.2251 0.2327 (31) 
Example 2 : [/}2°«'ds 0.0742 0.0730 


Example3 :f) 2 n(w + 1)ds 0.3055 0.3053 


4 Error analysis 


In this section, our analysis shows that when using CHFs to approximate a 
function, the order of accuracy is O(h*). Let us approximate a function u(t), 
as (10), where 


un(t)= So ulih)di(t), n=3K, KEN. (32) 
i=0 
In the first step, for t € (jh, (j + 1h), J = 0, 3, 6,...,n — 3, using (5)-(9) 
and doing some computation, we obtain 


Un(t) = Do u(jh)ei(t) 


° 


= 6, (u(ih) + djalOu(U 4D) 
+ dysa(t)u((j + 2)h) + dj+0(tu (Gi +3)h) 

= win) (CU DNE- GBM EW G+ 999) 
rugs ny (CAME G +I G+) 
pugs 2n) (CHINE GDA E= G+8)0) 
tu(ih ay (2=INE- G4 DMEM G+NY 


Then by simplifying the current relationship, we have 
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#8 3 _ 2 2 a 7 3 
vai =n ( jh)” — 6h(t ny 11h? (t — jh) a) 


(t — jh)? — 5A(t — jh)? + 6h? (t — jh) 
TT u(jh he h) ( Dh 
. (t — jh)® — h(t — jh)? + 3h? (t — jh) 
+u(jh + 2h) ( 58 
G=9hy Sh Gh) +o GS ph) 
+ u(jh + 3h) ( 6h3 : 


Therefore 


Un(t) = uGh) 


5h) ez + 18u(jh + h) le + 2h) + 2u(jh+ aH) 
(agi (2 5u(jh +h) + 4u(jh + 2h) — u(jh4 n)) 

p) h? 
(e=9n)° (= + 8u(jh +h) — 3u(jh + 2h) + u(jh + 0) 

6 h3 


It is known that the kth, k = 1,2,3, order derivative of u(t) about the point 
(jh) is as follows: 


—1lu(jh) + 18u(jh + h) — 9u(jh + 2h) + 2u(jh + 3h) 


ul(jh) = an ee 
Pepe Qu(jh) — Bu(jh +h) + sulsh + 2h) — u(jh + 3h) + 0(h4, 
ot i) A NGS SG ABET) os) 


So, assuming h > 0 results in 


t— jh)” t—jh)° 
un(E) = u(jh) + (t= shu (Gh) + 4 2 Yuin) +! 2 Y u"(jh) +O(H4). 
(33) 
Expanding u(t) in the Taylor’s series, about the point t = jh, we have 


3 -7\k 
u(t) = 5- Cay WW) (jh) + O(t — jh)}. (34) 
k=0 , 


According to (33) and (34), we can obtain the error between the exact and 
approximate values of u(t) as follows: 
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u(t) — un(t) = O(t — jh)*. (35) 


Thus, for t > jh, 7 = 0, 3, 6,...,2—3 and h > 0, since (t — jh) < h, from 


(35), we attain 
u(t) — un(t)| = O(h*). (36) 


In the second step, for 7 = 1, 4, 7,...,n-2and jh <t< (j+J)h, we get 


j—1(t)u ((j — Ih) + G3 (thu Gh) 
+ bjsi(t)u((9 + 1)A) + Oj +2(t)u ((9 + 2)h) 


l| 
oo 


= ogh= my (E=E- Ge DME G+) 
ucin, (2G DNEW G4 DAE G+ 9)9) 
ruin ny (2=G= DME HVE G4 999) 
pug gany (C2 G= DAE ME GDH) 


Hence 


Un(t) = ugh) 


L(t jh) (= —h) — 3u(jh) + Gu(jh + h) — u(jh+ =) 


6h 
, fe jh)” (u(jh —h) — Qu(jh) + u(jh +h) 
2 h? 
tt jh)® (—u(jh — h) + 3u(jh) — 3u(jh + h) + u(jh + 2h) 
6 h3 
(37) 
As a reminder, the derivatives of u(t) about the point jh are as follows: 
—2u(jh — h) — 3u(j jh +h) —u(jh + 2h 
jh —h) — 2u(jh h+h 
1! Gn) = MIR=M) = BUG) FMAM | oa) 
—u(jh — jh) — Ih+h jh + 2h 


As h — 0, from (37)-(38), we get 
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i= 3h) t—jh)° 
uin(t) = u(jh) +(t—gh)ul (jh) + 4 2 Y utt(gn) +! : Y u"(jh) +O(H4). 
(39) 
Considering (39), (34), and h > 0 for j = 1, 4, 7,. — 2, one has 


|u(t) — un(t)| = O(n"). (40) 


In the final step, for 7 = 2, 5, 8,...,n -landte (jh, (7 + 1)h), we attain 


= Sl uGihoilt 
i=0 


= $j-2(t)u ((j — 2)h) + 4)-1(8u (i Dh) 
+ &y(®u ik) + dja (8) (i+ DP) 
= vin 29) (E= GDA ERE GDN) 


—6h3 


Lith m (© G me Ie aay 


et é = G= 2) t= GDM = G+ yi) 


(t(j —2)h) (tJ —1)h) (tn) 
( i ) 


As a result, 


un(t) = u(jh) + (t — jh) (““ 2h) — 6u(gh = Sulsh) + 2u(Jh +1) 
(t— jh)” (ut —h) — 2u(jh) + u(jh + ») 

2 h2 

ie Ga aoc nae ee 4) 

6 h3 ; 


(41) 


On the other hand, according to the derivatives of u(t) about the point jh, 
(41) can be written as follows: 


Ee 


7 3 
2 ) ul" (jh) O(h*). 
(42) 

Thus, assuming (42), (34), and h > 0, for 7 = 2, 5, 8,...,n —1, one has 


uin(t) = ugh) + (t— gh)ul(Gh) + IN wrgny + § 


|u(t) — un(t)| = O(h*). (43) 
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Finally, for t € (jh, (g + Dh), 7 = 0,1,2,...,n, and h > 0, using (36), (40), 
and (43), we get 
|u(t) — un(t)| = O(n’). (44) 


5 Numerical algorithm 


In this section, a numerical algorithm is offered to solve problem (1). Consider 
the following nonlinear fractional integral-differential equation with weakly 
singular kernel : 


© u(t) = g(t) + p(t)u(t) +r ‘J (t — s)Pu™(s)ds, (45) 
a>0, 0<8<1, meEN, teETI(t). 


First, putting —8 =w—1,0<w <1 in the third term on the right of this 
equation, we get 


E (ts) Pu™(s)ds = T(w) (a if (t — s)*"w"(a)ds) , S2e<1 


(w 


By the definition of Riemann-—Liouville fractional integral operator, [24], the 
current relationship can be rewritten as follows: 


| (t — s) Pu (s)ds = T(w) I? (u™(t)) . (47) 


Now, by applying (47), the Riemann—Liouville integral operator of order a 
on the both sides of (45), one gets 


Ip (9 Deu(t)) = I? (g(t) + LF (w()u(t)) + ADwYETP (u™ (4) , 
u(t) = z(t) + IP (g(t) + IP (p@)u(t)) + AT (we (u(t), (48) 


where 


fa]—1 
z(t) = » u (0) — 


Now, by approximating the functions in (48) by CHFs (10) and (19), we 
attain 


u(t) ~ S~ aigi(t) = AT®(t), Am = [a0,41,---@n]”. (49) 
1=0 

u™(t) = > ay b(t) = Al &(t), Am = [ag’, aT’,... jal (50) 
1=0 
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a(t) ~ > 2(th)gi(t) = 27 H(t), Z = [2(0),2(h),-..,2(n)]7, (51) 


g(t) = Dd 9ih) ill) =GTW(t), G=[9(0),g(h),.-.,g(nh)]", (52) 


p(t) > Plih)di(t) =p’ W(t), P= [p(0),p(h),...,p(nh)]”, (53) 


wherein n is an integer multiple of 3. Utilizing (20)-(22), and (18) and the 
substitution of (49)—(53) in (48) become 


ATO(t) = ZT H(t) + 12 (GT O(t)) + [2 (PT S(t) @(t)7 A) + AT (W/E TY (AT G()), 
ATO(t) — Z7 &(t) — IP (GT G(t)) — If (PT diag (®(t)) A) — Ww) let’ (AF ®(t)) = 0, 
AT O(t) — ZT @(t) — 12 (GT G(t)) — [2 (P7 diag (A) ®(t)) — AT(w) IS + (AT &(t)) = 0, 
AT O(t) — Z7 @(t) — GTQMa(t) — PT diag (A) Q( &(t) — AT(w) AT Qt) B(£) = 0. 
Thus 


A? — ZF — G™Q@ — PT diag (A) Q™ — AT (w) AZ. Qt”) = 0, 
a>0, O0<w<1, w=1-8. (54) 


This system has the dimension (n+ 1) x (n+ 1). 
Suppose that 


Qe — eae Q’ = (9: ’ 1,9 =0,1,2,...,n. (55) 
Then, from the operational matrix (23), one gets 


Ye =H =O, 120, 1 2jacn gts 
Vij = ij = 0, j =1,3,...,n—-1, @#=794+3,j9+4,...,n, 
Vj = 915 = 0, j = 2,4,...,7, 4=94+2,74+3,...,n. 


Using (54), the unknown coefficients can be determined. We start to find the 
first unknown coefficient as follows: 


ay = z(0). (56) 


In the next step, we get 
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aa ~ [2(14)] ~ ] 35 oh) | 


i=0 


| 
oe (th) ya; | — |AT(w) % bia" | = 0, 
soon | tal = I~ atin) a | 
os p(ih) rans ~ | AP Ww) ¥ baci” =0, 
fag] ~ fe) ~ | oth) | 
z “fe plihyiaa | : far) 3 a0" = 


Solving the first system allows us to calculate a , a2, and a3, then we solve 


the following system: 


|— Fe(4ny] — | oti | 


(ih) isa | fares acs" =i, 
( 


Vw 


= [e(6h)] = [52 arr] 


eqs : | 

7 | 0 
Sshewiet eqs : [as 2 
y 2 D> p ih) 7 [pre yy 60s" 256) 
cas: [as] ~ [(6%)] ~ [32 oth) 


1=0 


= [= p(ih) soa | = [arw) p> 6a" =0. 


By solving system 2, the values of the unknown parameters a4, a5, and ag 
are calculated. Then we can get the values of a7, ag, and ag using system 


3: 


a6 


car: or] ~ 2(7h)] ~ [3° oi re] 
: 
do Plih) vrai yu Dra,” | = 0, 
E et =| eo ih . | 
ee) le ek — [ariey 3 dae" | =0 
can: [an] ~ [2(90)] — [5- ary 
_ [3 plan) io | _ [arw) » 2 ia" = 0. 
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The process can be continued up to the following form: 


vi g(th) Vi(n—2) 


etna [ana] ~ (ln 2)8)] — | 


_ bs 0-9 | - [arw) Y) Pin—ayas™ | = 0, 
b i=0 


edn-1 [en-1) ~[e((n-— Dh) LE slit) 
= = lib) on—ye | = [arw) ps Ayre =f 
cam an] ~ 2(nh)] — [3° othr 


i=0 


2 D> plihyrinas 7 [arw) » 8inas™ =9) 


system n/3: 


As a result, the values of an—2, Gn—1, and a, are derived using system n/3. 
Therefore, we can obtain an approximate solution via (10). To solve the 
nonlinear equations, see [34]. The computations were handled by MATLAB 
package. The following theorem outlines the proposed method. 


Theorem 2. Consider the principal problem (1). To obtain a numerical 
solution to (1) using CHFs, the following iterative algorithm is offered: 


Proof. See the scheme proposed in this section, (56)—(57). 


6 Convergence analysis 


In this section, we will verify the convergence of the numerical proposed 
scheme. 


Theorem 3. Let u,(t) be the numerical solution of (1) obtained by the 
proposed method in Section (5). Moreover, u(t) is an exact solution and 
E,,(t) is the residual error for numerical solution. Also, suppose that M and 
K are positive constants. Then, E,,(¢) tends to zero, as n — oo, where 


M == sup Ir-"(ay(t = 7)*'p(7) 
t,r€[0,T] 


K= = sup mL Ma +w)(t— san 
t,r€ [0,7] 


: | 


Proof. Applying the Riemann-—Liouville integral operator of order a and (48), 
it is appropriate to rewrite (1) in the integral form 


u(t) = 2(t) + IP (g(t) + IP (p@ult)) + AT(w) et? (w(t), (57) 


where 
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Step 1: Inputs, n (integer multiple of 3), a, 6, A,T, g(t), p(t), 
w™(0), i=0,1,...,fa]—1. 
Step 2: Setw=1-—8, h=T/n, and t; =ih, i=0,...,n 


n—-1 ; 
Step 3: z(t) = 0 uO (0)4. 


1=0 
Step 4: Compute the elements of Q( = [yj] and Q(t”) = [6;)], 
i,j =0,..., 
Step 5: Set and solve recursive trivariable system v, v = 1: n/3. 
ag = 2(0), 
for v=1:n/3 


Solution of the vth system determines 
the unknown parameter. 


[cao] — f((80 —2)0)] = | oih)ruea0-a | 
-|= p(ih)¥i(3v—2)4 a 
| P) ¥ Bcan-2yae” 
[ace] ~ f((80 ~ 19h) = [5° 9th) r4e0-1 | 
system: ) — [3° p(ih) rans 
- fre) 3 csr” | = 0, 
[any] ~ Fe(30)A)] ~ [5° ality 
: > lib 
Ss 
- fre) 3 cues” = 0, 


= 0, 


_— | 


end. 
Step 6: Calculate fully a;, 7=0,1,...,n 


Step 7: Define CHFs: (¢;(t), «=0,1,...,n). 


n 
Step 8: Determine the approximate solutions: un(t) = >> aidi(t). 
i=0 


Thus, u,(t) satisfies the following equation: 


—, a>0, w=1-6, 0<8<1, teEl(t). 
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Un(t) = z(t) + Te (9(t)) + Le (v(t)un(t)) + AT W)ET? (u(t). (58) 


If the residual function E,,(t) is not zero, then we can obtain it by using the 
following relation: 


En(t) = enlul(t) — Jelul(t) — Ver] (0), (59) 
where 
enlul(t) = ult) ~ unl), (60) 
1 : a-—l 
Fel) = Fay f - 1) wa Yulr) —tinr))ar, (61) 
and 
vertu) = ey f wm) —uer)ar. (62) 
Then, we get 
En (| < lenled(t)l + ated el + vet fer] (a). (63) 


For t € (ih, («+ 3)h),i = 0, 3, 6,...,n — 3, according to (44), the approxi- 
mation of the absolute error using CHFs yields 


|u(t) — un(t)| = O(n"). (64) 


By using (60), we have 
len [u](t)| = O(h*), (65) 


when h > 0, |en[u](t)| + 0. Then, by using (61) and (64), we attain 


FeO = roy ff = Pe) ule) — tar) 
< of (¢—7)°* |u(r) — ua(r)| ar 
< ee [b= 7)" u(r) = un(r)| ar 
< MO(h*), (66) 


wherein 


M= sup |I7l(a)(t— 7) p(n) and whenh —> 0, | J2-[ul(t)| > 0. 
t,r€[0,T] 


In addition, the following inequality holds [11]: 
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Ju (t) — un (t)| < mL u(t) — un(é)|, (67) 


where L = |(max(u(t), wn (t)))"]. As well, from (62), (64), and (67), we 
have 


vel = aoa f e-em) - una 
|A| ‘ atw—1 m m 
<aeay ff Oo) — ear 
mL |A| : a+tw-1 
Fray f, C—O ue) ule er 
< KO(h*), (68) 


wherein 


K= sup |AmLT~'(a+w)(t— a | andash—>0, |V,0**[u'™](t)| > 0. 
t,rE[0,T] 


Then, from relations (65), (66), (68), and (63), it is obvious that |E,,(t)| tends 
to zero, as h + 0, or n > co. 


7 Numerical examples 


In this section, the theoretical results of the previous sections are used for 
solving linear and nonlinear fractional integral-differential equations with the 
weakly singular kernel, that is, the initial condition equation (1). For assess- 
ing the accuracy of the scheme, let us define the maximum absolute error 
(LZoo-norm error) as 


Intl = sup {|u(ti) — un(te) |}. (69) 


‘ts =ih] e=0 


Using this definition, the order of convergence, with respect to this norm, is 
introduced as follows: 


Rate = log, ( ee ) ’ (70) 


For some problems, there are no exact solutions, so the L2-norm error is 
calculated by the following formula: 


iS 


lénllo = (Se toute wont , t;=th, i=0,...,n, (71) 


i=0 
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where un(t), n = T/h is the approximate solution defined as (49). In addi- 
tion, the results of different values of @ are compared with each other and 
with MHFs method, [22]. 


Example 1. Consider a nonlinear integral-differential equation with weakly 
singular kernel [22]: 


CDPult) = alt) + 2(u(t) + | (t—s) °u?(s)ds, t € [0,1], 
2 Vat (7) \ 23 
t) = 3t 2 fF, t)=0, u(0)=0 
g(t) r(3) p(t) (0) 


For a = 1 and B = $, the exact solution is u(t) = t?. Approximate 
numerical results using different values of n are shown in Tables 1-3 and Fig- 
ures 2-4. Table 1 shows the approximate and exact solutions to the problem 
at some points. Also, the L,.-norm errors and convergence orders obtained 
by the current method are compared with the MHFs method [22] and the 
methods presented in [14, 5] in Table 2. Table 3 shows the comparison of 


the result of the /g-norm error ||&,||, obtained by the proposed method. It is 
clear from Table 2 that the results of the present method for less or similar n 
are better than the results obtained in [22]. Figure 2 indicates the behavior 
of absolute errors for Example 1. Also, Figure 3 shows the logarithm of the 
L.-norm errors. As can be seen from the plot, as n increases, the error 
decreases. Also, the comparison of the results obtained for different values 
of alpha with the exact solutions of the equation is plotted in Figure 4. 


Table 1: Numerical results of Example 1 
Points Exact Approximate Approximate Approximate 
solutions solutions solutions solutions 

8 u(s) n= 12 n = 24 n = 48 
0 0.00000000000 | 0.00000000000 | 0.00000000000 | 0.00000000000 
Tho 0.00057870370 | 0.00057871890 | 0.00057870463 | 0.00057870368 
M% 0.00462962963 | 0.00462979683 | 0.00462962554 | 0.00462962982 
Ty, 0.01562500000 | 0.01562510533 | 0.01562500074 | 0.01562499998 
Ti, 0.03703703704 | 0.03703629382 | 0.03703707090 | 0.03703703603 
She 0.07233796296 | 0.07233970584 | 0.07233791123 | 0.07233796502 
Ty 0.12500000000 | 0.12500014128 | 0.12499999725 | 0.12499999984 
Th 0.19849537037 | 0.19849178055 | 0.19849551759 | 0.19849536600 
2) 0.29629629630 | 0.29630236028 | 0.29629611403 | 0.29629630318 
3, 0.42187500000 | 0.42187499884 | 0.42187498758 | 0.42187499953 
% 0.57870370370 | 0.57869407976 | 0.57870407228 | 0.57870369262 
Thy 0.77025462963 | 0.77026833822 | 0.77025420450 | 0.77025464509 
iL. 1.00000000000 | 0.99999965644 | 0.99999997045 | 0.99999999902 
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Table 2: Comparison of the Loo-norm error and convergence order for Example 1 


MHFs method [22| CHFs method 
n En lle Rate of convergence |_ n lEnlloo Rate of convergence 
4 | 4.81704E — 03 3.93 3 5.643508 — 03 4.22 
8 3.15569E — 04 3.49 6 3.03330E — 04 4.47 
16 | 2.81379E — 05 3.88 12 1.37086. — 05 4.71 
32 | 1.90987E — 06 3.93 24 | 5.24867 — 07 4.85 
64 | 1.25460E — 07 3.95 48 1.82605E — 08 4.92 
128 | 8.09506 — 09 3.97 96 6.04389 — 10 4.96 
256 | 5.16785E — 10 = 192 | 1.950051F — 11 = 


Method presented in [14] | Method presented in [5] 
4 3.52F — 09 3.58 — 04 
8 1.40£ — 14 1.11022F — 16 
16 1.51f — 14 1.11022F — 16 


Table 3: Numerical results of the L2-norm error functions ilEnlly for Example 1 


n 3 6 12 24 48 96 192 
lEnllp | 7.9E — 03 | 3.5E—04 | 1.9F—05 | 9.3E—07 | 44E— 08 | 2.0F—09 | 96-11 


x10° 


Error for m=12 \ 
Error for m=24 x10" \ 
1.2 |—¥— Error for m=48 


Absolute error 
eo = my wb 
eS 


0 005 01 015 02 | 


Absolute error 


Figure 2: Absolute errors of Example 1, for n = 12, 24,48 


Example 2. Consider the following nonlinear fractional integral-differential 
equation with weakly singular kernel [22] : 


6 Dru(t) = g(t) + 
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-norm error) 


90 


log, (| 


-35 


Figure 3: Logarithm of the Loo-norm error in Example 1 


Exact solutions a=1 

— ~~ Approximate solutions a=1 
| |———Approximate solutions «=0.95 
— ~~ Approximate solutions a=0.85 
1.2 | |— — ~ Approximate solutions a=0.75 ie 7 


0 01 


Figure 4: Exact and approximate solutions of Example 1, for n = 12 


For a = § and 6 = 35, the analytic solution is u(t) = t?. Approximate 
numerical results using different values of n are shown in Tables 4-6, and Fig- 
ures 5-7. Table 4 shows the approximate and exact solutions to the problem 
at some points. Table 5 indicates the L,,-norm errors and convergence orders 
for various values of n. Table 6 shows the L2-norm errors at some values of n. 
As can be compared in Table 5, this new method provides a higher order of 
convergence compared to the other method. Figure 5 indicates the absolute 
errors of Example 2, at n = 12,24,48. Figure 6 shows that the logarithm 
of the L,,-norm error decreases as n increases. A comparison between the 
changes in the fractional orders of the equation is shown in Figure 7. 


Example 3. Consider the following linear fractional integral-differential 
equation with weakly singular kernel [22]: 


CS D@u(t) = g(t) + p(t)u(t) +| (t— s) Pu(s)ds, t € (0, 1], 
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Table 4: Numerical results of Example 2 


Points Exact Approximate | Approximate | Approximate 
solutions solutions solutions solutions 
8 u(s) n=12 n = 24 n= 48 

0 0.00000000 0.00000000 0.00000000 0.00000000 
Th 0.02405626 0.02340582 0.02394753 0.02400902 
TM 0.06804138 0.06771773 0.06790284 0.06800459 
Ty, 0.12500000 0.12449685 0.12487612 0.12496621 
Ti, 0.19245009 0.19199415 0.19232923 0.19241677 
whe 0.26895718 0.26849729 0.26883160 0.26892254 
Ty 0.35355339 0.35305089 0.35341692 0.35351566 
The 0.44552819 0.44495789 0.44537319 0.44548528 
/3 0.54433105 0.54365866 0.54414702 0.54428012 
3), 0.64951905 0.64868490 0.64929131 0.64945599 
He 0.76072577 0.75964752 0.76043179 0.76064433 
lh 0.87764152 0.87619255 0.87724539 0.87753181 
iL 1.00000000 0.99796079 0.99944320 0.99984577 


Table 5: Comparison of the Loo-norm error and convergence order for Example 2 
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MHFs method [22] CHFs method 
n lEnlloo Rate of convergence | n En lloo Rate of convergence 
4 | 1.39991F — 02 2.08 3 | 4.05382E — 02 2.38 
8 | 3.30803E — 03 1.90 6 | 7.80949E — 03 1.94 
16 | 8.84780E — 04 1.84 12 | 2.03921E — 03 1.87 
32 | 2.46509F — 04 1.83 24 | 5.56795E — 04 1.85 
64 | 6.92324E — 05 1.51 48 | 1.54232F — 04 1.84 
128 | 2.42498F — 05 1.50 96 | 4.30102E — 05 1.84 
256 | 8.57216E — 06 - 192 | 1.20325E — 05 - 


Table 6: Numerical results of the L2-norm error functions ilénlly for Example 2 


n 3 6 12 24 48 96 192 
lEnllp | 3-56-02 | 7.2E—03 | 2.4E—03 | 856 —04 | 33E—04 | 13E—04 | 5.0F—0 
32 1 
t)=-—t?2, u(0)=0. 
rt)=—st#, u(0) 


For a = } and = 3, the analytic solution is u(t) = t? + t3. Tables 7-9 
and Figures 8-10 show approximate numerical results using different values 
of n. Table 7 indicates the approximate and exact solutions to the problem 
at some grid points. Table 7 shows the advantage of the proposed method 
compared to the MHF method by presenting the order of convergence and 
the maximum norm error. Figure 8 shows the behavior of absolute errors 
for Example 3. Figure 9 shows the logarithm of the Z.-norm errors. In 
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2.5 T , r r , r T T T 

Error for m=12 
Error for m=24 
——*— Error for m=48 


nD 


Absolute error 


15 


Absolute error 


log, (I, -norm error) 
Ss 


Figure 6: Logarithm of the L..-norm error in Example 2 


45 r r r r r r : : : 
Exact solutions a=2/3 
4+ |— — — Approximate solutions «=2/3 1 
— ~~ Approximate solutions o=3/4 f 
3.5 | |— — — Approximate solutions o=7/12 4 
— —~Approximate solutions a=1/2 ; 
3h | 
/ 
25F 4 
_ / 
ol ge | 


Figure 7: Exact and approximate solutions of Example 2,for n = 12 
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addition, a comparison of the results for different values of a with the exact 
solution of the equation is shown in Figure 10. 


Table 7: Numerical results of Example 3 


Points Exact Approximate | Approximate | Approximate 
solutions solutions solutions solutions 
8 u(s) n=12 n = 24 n = 48 

0 0.00000000 0.00000000 0.00000000 0.00000000 
Tho 0.03697789 0.03708156 0.03699820 0.03698362 
% 0.09634983 0.09644877 0.09637681 0.09635574 
va 0.17311513 0.17323059 0.17314347 0.17312128 
13 0.26815746 0.26828582 0.26818799 0.26816396 
He 0.38354663 0.38368969 0.38357860 0.38355352 
Th 0.52185026 0.52199562 0.52188416 0.52185755 
he 0.68589934 0.68605454 0.68593554 0.68590705 
2) 0.87868327 0.87885234 0.87872136 0.87869145 
3, 1.10329522 1.10346969 1.10333552 1.10330386 
% 1.36290039 1.36308529 1.36294318 1.36290952 
Th 1.66071634 1.66091514 1.66076132 1.66072598 
1 2.00000000 2.00020619 2.00004743 2.00001015 


Table 8: Comparison of the Loo-norm error and convergence order for Example 3 


MHFs method [22] CHFs method 
n lEnlloo Rate of convergence | n En lloo Rate of convergence 
4 | 1.14967E — 03 1.86 3 | 1.68472E — 03 1.20 
8 | 3.15569E — 04 1.83 6 | 7.35848E — 04 1.84 
16 | 8.90424EF — 05 2.20 12 | 2.06192F — 04 212. 
32 | 1.93665E — 05 2.30 24 | 4.74299F — 05 2.22 
64 | 3.94225E — 06 2.09 48 | 1.01541E — 05 2.27 
128 | 9.26153E — 07 2.13 96 | 2.10201E — 06 2.30 
256 | 2.12088E — 07 - 192 | 4.27476E — 07 - 
Table 9: Numerical results of the L2-norm error ||En||, for Example 3 
n 3 6 12 24 48 96 192 
lfnllp_ | 12H-03 | 10E-—03 | 4.2E—04 | 14-04 | 42-05 | 1.2L—05 | 3.6E —06 


8 Conclusion 


In this paper, we proposed a numerical scheme for solving a class of non- 
linear fractional integral-differential equations with weakly singular kernels 
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Exact solutions a=1/3 
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Figure 10: Exact and approximate solutions of Example 3, for n = 12 


based on CHF'’s. CHF and the corresponding operational matrix were in- 
troduced. The proposed method transforms the original problem into an 
iterative algorithm, including polynomial equations with three unknown co- 
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efficients, using the fractional-order operational matrix of integration. An 
analysis of the method’s absolute errors and convergence was conducted. In 
order to validate the accuracy and effectiveness of this new method, three 
numerical examples were presented. In Example 1, the absolute error is lower 
at the nodal points near the beginning of the interval, as shown in Figure 
2. Table 2 shows that the new method offers more accurate solutions at the 
same lengths h than the MHFs approaches. In Examples 2 and 3, the error 
clearly increases as the time variable approaches one; see Figures 5 and 8, 
respectively. A study of the results shows that, generally, as n increases, 
the accuracy of the approximate solution increases, and the absolute error 
decreases. One of the advantages of this proposed algorithm is that instead 
of solving a system of (n+ 1) x (n+ 1) equations, it only needs to solve n/3 
systems of three-variable nonlinear equations. In addition, the order of con- 
vergence for the cubic hat functions is O(h*), while the order of convergence 
for the generalized hat functions method [6] and the MHFs method [22] are 
O(h?) and O(h3), respectively. Finally, the proposed method (CHFs) can be 
used for a large number of similar problems, and we will continue to work on 
developing this method. 
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